Objetivo

Realizar mapas con la mayor granularidad posible, a escala intraurbana, de vulnerabilidad social ante la pandemia de COVID-19.

Con ello se busca contribuir en la toma de decisiones por parte de entidades municipales para la planificación y distribución de atención a nivel intraurbano (víveres con bancos de comida, atención médica in situ, etc.).

En esta instancia se utilizarán solo fuentes abiertas con cvertura global, de modo que el procedimeinto sea reproducible en cualquer ciudad de Latinoamérica sin requerir fuentes locales específicas. Por supuesto, los datos locales de gran valor y podrán ser incorporados a éste procedimiento para refinar o extender los resultados.

Esta prueba de concepto toma a Quito, Ecuador, como ubicación a analizar.

Fuentes de datos

Paquetes a utilizar

library(tidyverse) # funciones útiles en general para manipulación y visualización
library(sf) # para datos georeferenciados en formato vectorial
library(osmdata) # para acceso a datos de OpenStreetMap
library(leaflet) # para visualización de datos geográficos
library(nngeo) # para identificar elementos próximos en el espacio  
               # No disponible en CRAN, usar devtools::install_github("michaeldorman/nngeo")

Descarga de datasets

Estimados demográficos

Población general

Descargamos los datos de estimados demográficos para Ecuador producidos por Facebook

url <- "https://data.humdata.org/dataset/58c3ac3f-febd-4222-8969-59c0fe0e7a0d/resource/c05a3c81-a78c-4e6c-ac05-de1316d4ba12/download/population_ecu_2018-10-01.csv.zip"

zipfile <- tempfile()
download.file(url, zipfile)
filename <- unzip(zipfile, list = TRUE)$Name

densidad <- read_csv(unz(zipfile, filename))

unlink(temp)

El dataset contiene puntos definidos por pares de coordenadas en proyección Mercator, y estimados de población en torno a esa posición en 2015 y 2020, para todo Ecuador:

head(poblacion)

Cada par latitud/longitud representa un área que corresponde a 1 segundo de arco de resolución (aproximadamente 30m x 30m)

Población mayor a 60 años

De la misma fuente descargamos la cantidad estimada de personas de más de 60 años

url <- "https://data.humdata.org/dataset/58c3ac3f-febd-4222-8969-59c0fe0e7a0d/resource/904d8988-d18d-41a5-a7f7-22668204cefe/download/ecu_elderly_60_plus_2019-06-01_csv.zip"

zipfile <- tempfile()
download.file(url, zipfile)
filename <- unzip(zipfile, list = TRUE)$Name

personas_mayores <- read_csv(unz(zipfile, filename))

unlink(temp)

En este caso el dataset contiene el estimado de personas mayores a 60 años viviendo en el área que corresponde a cada punto:

head(personas_mayores)

Límites de la ciudad

El polígono de las fronteras administrativas de la ciudad puede encontrarse en Nominatim: https://nominatim.openstreetmap.org/

Debido a la infinidad de marcos territoriales y legales que existen en el mundo para definir las fronteras de las ciudades, en general existen varias formas de interpretar cuales son los límites de una ciudad.

Nominatim es una base de datos global de nombres propios de lugares. Al realizarse una búsqueda, por ejemplo por “Quito”, se encuentran distintas entidades geográficas con ese nombre. La interfaz de Nominatim permite comparar opciones y verificar los límites en el mapa, identificando la entidad buscada.

En este caso, buscando por “Quito, Ecuador”, el resultado deseado es el quinto resultado: Quito, Pichincha, Ecuador (County), que coincide con los límites visibles en el SIG de la secretaría de Territorio de QUito

knitr::include_graphics("img/secretaria_territorio_quito.png")

Identificado el nombre, descargamos los límites.

ciudad <- "Quito, Ecuador"
# La posición en que aparece entre los resultados de Nominatim 
posicion <- 5

url <- URLencode(paste0("https://nominatim.openstreetmap.org/search.php?q=",
                        ciudad, "&polygon_geojson=1&format=geojson"))

destfile = paste0(tempdir(), "/limits.json")

download.file(url, destfile)

limites <- st_read(destfile)[posicion,]
## Reading layer `limitesquito' from data source `/home/havb/Downloads/limitesquito.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -78.94768 ymin: -0.5956779 xmax: -78.1649 ymax: 0.2562735
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
leaflet(limites) %>%
    addProviderTiles(providers$CartoDB.Positron) %>% 
    addPolygons()

Puntos de Interés

Para obtener datos georeferenciados locales realizaremos consultas a Overpass (http://overpass-api.de/), una interfaz que permite extraer información de la base de datos global de OpenStreetMap.

Overpass requiere que se especifique una “bounding box”, es decir las coordenadas de un rectángulo que abarque la zona de interés (en la práctica, los valores máximos y mínimos de latitud y longitud).

La generamos:

bbox <- matrix(st_bbox(limites), nrow = 2, dimnames = list(c("x", "y"), c("min", "max")))

bbox
##           min         max
## x -78.9476819 -78.1649035
## y  -0.5956779   0.2562735

Para descargar entidades georeferenciadas se requiere conocer las palabras clave con las que se identifican los registros en la base de OSM. Existe gran detalle para el tipo de datos georeferenciados disponbles: áreas de parques públicos, posición de oficinas de correo o cajeros automáticos, vías de ferrocarril, etc. La nomenclatura se puede consultar en https://wiki.openstreetmap.org/wiki/Map_Features

En este caso vamos a solicitar todas las vías de circulación (calles, avenidas, autopistas, etc) de la ciudad. En la base de datos de OSM ciertas entidades agrupadas bajo la categoría “Healthcare” resultan de inmediato interés:

key value description
amenity clinic A medium-sized medical facility or health centre.
amenity hospital A hospital providing in-patient medical treatment
amenity pharmacy Pharmacy: a shop where a pharmacist sells medications

Para esta prueba de concepto descargaremos clínicas y hospitales.

clinicas <-  opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "clinic") %>% 
  osmdata_sf() 

hospitales <-  opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "hospital") %>% 
  osmdata_sf()

Dependiendo del nivel de detalle con el que hayan sido registrados, la mayoría de las clínicas y hospitales están representadas por puntos, pero en algunos casos por los polígonos de su planta. Para nuestros fines sólo necesitamos los puntos. Convertiremos los polígonos en puntos tomando su centroide.

todo_a_puntos <- function(osm_query_sf) {
  
  puntos <- select(osm_query_sf$osm_points, "name")
  centroides_poly <- st_centroid(select(osm_query_sf$osm_polygons, "name"))
  centroides_multipoly <- st_centroid(select(osm_query_sf$osm_multipolygons, "name"))

  todo <- rbind(puntos, centroides_poly, centroides_multipoly)
  
  # Solo retenemos los puntos con nombre, ya que los que tienen ese campo vacío
  # parecen ser ubicaciones repetidas
  filter(todo, !is.na(name))
}
  
hospitales <- cbind(todo_a_puntos(hospitales), clase = "hospital")
clinicas <- cbind(todo_a_puntos(clinicas), clase = "clinica")

hospitales_y_clinicas <- rbind(clinicas, hospitales)
## Reading layer `hospitales_y_clinicas_quito' from data source `/home/havb/Downloads/hospitales_y_clinicas_quito.geojson' using driver `GeoJSON'
## Simple feature collection with 457 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -78.89925 ymin: -0.5168509 xmax: -78.20187 ymax: 0.2484693
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

El resultado:

hospitales_y_clinicas
leaflet(hospitales_y_clinicas) %>%
    addProviderTiles(providers$OpenStreetMap) %>% 
    addMarkers(popup = ~paste(clase, name)) %>% 
    addPolygons(data = limites, fill = NA)

Procesamiento

División de la ciudad en zonas de análisis

Aquí es donde cada ciudad puede decidir su unidad geográfica de análisis. Lo ideal serían unidades estadísticas censales, con la mayor resolución (lo más pequeñas) que se disponga.

EL análisis a nivel de agregación similar al barrio (o distritos, comunas, etc) no es deseable debido a su baja resolución espacial. En ausencia de cartografía censal de alta granularidad, puede realizarse la partición de la superficie de la ciudad en celdas arbitrarias, lo suficientemente pequeñas.

Para este ejemplo particionaremos la superficie de la ciudad en celdas de unos 500m de radio (aproximando una superfice de 0.785 km2)

# Usamos una re-proyección equiareal para mayor precisión al calcular superficies

n_celdas <- limites %>% 
  st_transform(crs = "+proj=laea +ellps=WGS84 +units=m +no_defs") %>%
  st_area() %>% 
  {. / (pi * 500^2)} %>% 
  as.numeric() %>% 
  round()

n_celdas
## [1] 5364

Definiremos entonces 5364 celdas.

celdas <- limites %>%
  st_sample(size = 100000) %>%  # sampleamos al azar un número grande, a gusto
  st_coordinates() %>% # extraemos sendas columnas con long y lat
  kmeans(centers = n_celdas) %>%  # nuestros N grupos de puntos con separación máxima entre si
  .$centers %>% # sólo queremos los centroides
  as_tibble() %>% # convertimos en tibble que es lo que le gusta a sf
  st_as_sf(coords = c("X", "Y"), crs = 4326) %>% # pasamos los centroide a objeto espacial
  st_union() %>% # los combinamos en un sólo objeto multipunto
  st_voronoi() %>% # particionamos en voronoi el espacio donde estan
  st_collection_extract("POLYGON") %>% # del resultado extraemos los polígonos
  st_intersection(limites) %>% # recortamos el cuadrado de los voronoi de acuerdo a los límites
  st_sf %>% # Convertimos en dataframe espacial
  mutate(id = rev(row_number())) # agregamos columna con id
## Reading layer `celdas' from data source `/home/havb/Downloads/celdas.geojson' using driver `GeoJSON'
## Simple feature collection with 5364 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -78.94768 ymin: -0.5956779 xmax: -78.1649 ymax: 0.2562735
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ggplot(celdas) +
  geom_sf(aes(fill = NULL)) +
  theme_void()

Estimados demogŕaficos por área

Para el ejemplo usaremos los estimados de población mayor a 60 años, como grupo de riesgo. También podría considerase tanto la población general como los mayores, combinando ambas variables un un índice agregado.

En aras de la performance, primero extraemos del dataset nacional los datos que caen dentro de la bounding box de la ciudad (esta operación es muy rápida). Habiendo limitado así los puntos a clasificar, toma mucho menos tiempo verificar a que celda de la ciudad corresponde cada punto del dataset.

personas_mayores_quito <- personas_mayores %>% 
  filter(between(latitude, bbox["y", "min"], bbox["y", "max"]), 
         between(longitude, bbox["x", "min"], bbox["x", "max"])) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
  st_join(celdas, left = FALSE)
ggplot(personas_mayores_quito) +
  geom_sf(aes(color = population), size = .0001, alpha = .3) +
  scale_color_viridis_c() +
  theme_minimal() + 
  labs(title = "Población de mayores de 60 años por punto", color = "personas")

Agregamos la cantidad de personas por celda, en base a los punto que caen en cada una.

celdas <- personas_mayores_quito %>% 
  st_set_geometry(NULL) %>% 
  group_by(id) %>% 
  summarise(population = sum(population, na.rm = TRUE)) %>% 
  {left_join(celdas, .)}
## Reading layer `celdas_con_data' from data source `/home/havb/Downloads/celdas_con_data.geojson' using driver `GeoJSON'
## Simple feature collection with 5364 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -78.94768 ymin: -0.5956779 xmax: -78.1649 ymax: 0.2562735
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ggplot(celdas) +
  geom_sf(aes(fill = population), color = NA) +
  scale_fill_viridis_c() +
  theme_minimal() + 
  labs(title = "Población de mayores de 60 años por área", fill = "personas")

Efector de salud más cercano a cada área

Para estimar distancias entre cada celda y su efector de salud más cercano, crearemos una matriz Origen-Destino. Tomamos el centroide de cada celda, y determinamos el efector de salud (hospital o clínica) más cercano a su posición.

Como paso previo, descartamos celdas con población menor a un umbral dado:

umbral <- 1 # Solo celdas con al menos 1 persona

celdas <- celdas %>% 
  filter(population >= umbral)

Encontramos el elemento de la lista Y (sitios de salud) más cercano a cada elemento de la lista X (centroides de celdas).

En esta instancia tomaremos distancia lineal. En una futura iteración, estimaremos distancia a pie a través de la grilla de calles local usando OSRM

# Calculamos el centroide
centroides <- st_centroid(celdas)

# Identificamos el índice (la posición en el dataframe) de hospital/clínica más cercano a cada uno
id_cercanos <- unlist(st_nn(centroides, hospitales_y_clinicas))

# Calculamos la distancia entre cada par (por ahora lineal, a mejorar como distancia a pie por calles)
distancia_a_salud <- st_distance(centroides, hospitales_y_clinicas[id_cercanos,], by_element = TRUE)
# eliminamos la unidad (m)
distancia_a_salud <- as.numeric(distancia_a_salud)

# Agregamos as distancias
celdas <- cbind(celdas, distancia_a_salud) 
#st_write(celdas, "~/Downloads/celdas_con_data_distancia.geojson")
celdas <- st_read("~/Downloads/celdas_con_data_distancia.geojson")
## Reading layer `celdas_con_data_distancia' from data source `/home/havb/Downloads/celdas_con_data_distancia.geojson' using driver `GeoJSON'
## Simple feature collection with 2105 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -78.94076 ymin: -0.5017812 xmax: -78.22107 ymax: 0.2517594
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ggplot(celdas) +
  geom_sf(aes(fill = distancia_a_salud), color = NA) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() + 
  labs(title = "Distancia a efector de salud más cercano", fill = "m")

Identificacion de zonas de riesgo

Combinando los datos demográficos con los de cercanía a salud, resaltamos las zonas donde se concentra población mayor, y a la vez se verifican las mayores distancias.

ggplot(celdas) +
  geom_point(aes(x = population, y = distancia_a_salud)) +
  labs(x = "personas mayores a 60 años") +
  theme_minimal()

En base al gráfico de dispersión, consideraremos casos crítico a las celdas con más de 500 residentes mayores que distan más de 1.5 km del centro de salud más cercano

celdas <- celdas %>% 
  mutate(sitio_critico = population > 500 & distancia_a_salud > 1500)
ggplot(celdas) +
  geom_point(aes(x = population, y = distancia_a_salud)) +
  geom_point(data = filter(celdas, sitio_critico), 
             aes(x = population, y = distancia_a_salud), color = "red") +
  labs(x = "personas mayores a 60 años") +
  theme_minimal()

En el mapa:

celdas %>% 
  filter(sitio_critico) %>% 
  leaflet() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolygons(color = "red", popup = ~paste("población > 60 años:", round(population),"</br>",
                                            "distancia a salud:", round(distancia_a_salud), "m"))